Numerische Simulation Auf Massiv Parallelen Rechnern Wavelet Galerkin Schemes for 3d-bem Preprint-reihe Des Chemnitzer Sfb 393 Contents 1 Wavelets on Manifolds 2

This paper is intended to present wavelet Galerkin schemes for the boundary element method. Wavelet Galerkin schemes employ appropriate wavelet bases for the discretization of boundary integral operators. This yields quasisparse system matrices which can be compressed to O(NJ ) relevant matrix entries without compromising the accuracy of the underlying Galerkin scheme. Herein, O(NJ ) denotes the number of unknowns. The assembly of the compressed system matrix can be performed in O(NJ) operations. Therefore, we arrive at an algorithm which solves boundary integral equations within optimal complexity. By numerical experiments we provide results which corroborate the theory. Introduction Various problems in science and engineering can be formulated as boundary integral equations. In general, boundary integral equations are solved numerically by the boundary element method (BEM). For example, BEM is a favourable approach for the treatment of exterior boundary value problems. Nevertheless, traditional discretizations of integral equations su er from a major disadvantage. The associated system matrices are densely populated. Therefore, the complexity for solving such equations is at least O(N2 J), where NJ denotes the number of equations. This fact restricts the maximal size of the linear equations seriously. Modern methods for the fast solution of BEM reduce the complexity to a suboptimal rate, i.e., O(NJ log NJ), or even an optimal rate, i.e., O(NJ). Prominent examples for such methods are the fast multipole method [17], the panel clustering [20] or hierarchical matrices [19, 32]. As introduced by [1] and improved in [9, 10, 11, 12, 30], wavelet bases o er a further tool for the fast solution of integral equations. In fact, a Galerkin discretization with wavelet bases results in quasi-sparse matrices, i.e., the most matrix entries are negligible and can be treated as zero. Discarding these nonrelevant matrix entries is called matrix compression. It has been shown in [30] that only O(NJ) signi cant matrix entries remain. Concerning boundary integral equations, a strong e ort has been spent on the construction of appropriate wavelet bases on surfaces [7, 14, 15, 21, 26, 30]. In order to achieve the optimal complexity of the wavelet Galerkin scheme, wavelet bases are required with a su ciently large number of vanishing moments. Our realization is based on biorthogonal spline wavelets derived from the multiresolution developed in [4]. These wavelets are advantageous since the regularity of the duals is known [33]. Moreover, the duals are compactly supported which preserves the linear complexity of the fast wavelet transform also for its inverse. This is an important task for the coupling of FEM and BEM, cf. [22, 23]. Additionally, in view of the discretization of operators of positive order, for instance, the hypersingular operator, globally continuous wavelets are available [2, 5, 14, 21]. 1 The e cient computation of the relevant matrix coe cients turned out to be an important task for the successful application of the wavelet Galerkin method [13, 21, 27, 28, 30]. We present a fully discrete Galerkin scheme based on numerical quadrature. Supposing that the given manifold is piecewise analytic we can use a hp-quadrature scheme [21, 30] in combination with exponentially convergent quadrature rules. This yields an algorithm with asymptotically linear complexity without compromising the accuracy of the Galerkin scheme. The outline of the present paper is as follows. In section 1 we the construction of suitable wavelet bases on manifold. We treat the also globally continuous wavelet bases which are required for the discretization of boundary integral operators of positive order. With these bases at hand we are able to introduce the fully discrete wavelet Galerkin scheme in section 2. We survey on several practical aspects like setting up the compression pattern, assembling the system matrix and preconditioning. In section 3 we present numerical results which con rm the theoretical results quite well. The accuracy of the Galerkin scheme has never been deteriorated by the matrix compression. 1 Wavelets on manifolds For the treatment of boundary integral equations, wavelet bases de ned on manifolds are required which provide besides an approximation power also a certain number of vanishing moments. This by itself is a nontrivial task. We apply the idea of [14, 21]. For sake of brevity we recall only some basics of this construction. Moreover, we focus mainly on piecewise constant and linear functions. The outline is as follows. First, we introduce the general concept of a biorthogonal multiresolution analysis. Then, we de ne wavelets on the interval [0; 1] which are brought to the unit square by tensor product techniques. Utilizing a domain decomposition strategy, these wavelets are lifted onto the manifold via parametrization. 1.1 Wavelets and multiresolution analysis Multiresolution is by now a well-studied topic in signal processing. There are many excellent accounts about it, we refer the reader to the survey paper [6] and the references therein. Here we collect only some facts which are useful for our purpose. Let be a domain 2 Rn or manifold 2 Rn+1 . Then, in general, a multiresolution analysis consists of a nested family of nite dimensional subspaces Vj0 Vj0+1 Vj Vj+1 L2( ); (1) such that dimVj 2nj and [ j j0 Vj = L2( ): (2) 2 Each space Vj is de ned by a single scale basis j = f j;k : k 2 jg, i.e., Vj = spanf jg, where j denotes a suitable index set with cardinality j jj 2nj. A nal requirement is that these bases are uniformly stable, i.e., for any vector c 2 l2( j) holds kckl2( j) k jckL2( ) (3) uniformly in j. Furthermore, the single scale bases satisfy a locality condition diamsupp j;k 2 j: Instead of using only a single scale j the idea of wavelet concepts is to keep track to increment of information between two adjacent scales j and j + 1. Since Vj Vj+1 one decomposes Vj+1 = Vj Wj with some complementary space Wj, Wj \ Vj = f0g, not necessarily orthogonal to Vj. Of practical interest are the bases of the complementary spaces Wj in Vj+1 j = f j;k : k 2 rj = j+1 n jg: It is supposed that the collections j [ j are also uniformly stable bases of Vj+1. If = 1 [ j=j0 1 j; where j0 1 := j0, is a Riesz-basis of L2( ) we will call it a wavelet basis. We assume that these basis functions j;k are local with respect to the corresponding scale j, i.e., diam supp j;k 2 j and we will normalize them such k j;kkL2( ) 1. We note that at rst glance it would be very convenient to deal with a single orthonormal system of wavelets. But it was shown in [12, 30] that orthogonal wavelets are not completely appropriate for the e cient solution of boundary integral equations. For that reason we use biorthogonal wavelet bases. Then, we have also a biorthogonal, or dual, multiresolution analysis, i.e., dual single-scale bases e j = fe j;k : k 2 jg and wavelets e j = fe j;k : k 2 jg which are coupled to the primal ones via ( j; e j)L2( ) = I; ( j; e j)L2( ) = I: The associated spaces e Vj := spanfe jg and f Wj := spanfe jg satisfy Vj ? f Wj; e Vj ? Wj: (4) 3 Denoting likewise to the primal side e = 1 [ j=j0 1 e j; e j0 1 := e j0 ; then, every v 2 L2( ) has a representation v = 1 X j j0 1 X k2rj(v; j;k)L2( ) e j;k = 1 X j j0 1 X k2rj(v; e j;k)L2( ) j;k and moreover kvk2L2( ) 1 X j j0 1 X k2rj j(v; j;k)L2( )j2 1 X j j0 1 X k2rj j(v; e j;k)L2( )j2: We refer to [6] for further details. For the construction of multiresolution bases on domains and manifolds one is interested in the lter coe cients or mask coe cients associated with the scaling functions and the wavelets. Since additional boundary functions are introduced, these lter coe cients are not xed like in the stationary case. Therefore, we compute the full two scale relations j = j+1Mj;0; j = j+1Mj;1; and likewise for the dual counterparts e j = e j+1f Mj;0; e j = e j+1f Mj;1: We mention that these matrices will be banded and only the lter coe cients for some speci c scaling functions and wavelets have to be modi ed. That way, the advantages of the stationary and shift-invariant case are preserved as far as possible. 1.2 Biorthogonal spline multiresolution on the interval Our approach is based on the biorthogonal spline multiresolution on R developed by A. Cohen, I. Daubechies and J.-C. Feauveau [4]. These functions have several properties which make them favourite candidates for a wavelet Galerkin scheme discretizing boundary integral equations. The primal multiresolution consists of cardinal B-splines of the order d as scaling functions. Therefore, we have to deal with piecewise polynomials. This simpli es the construction of wavelet bases on manifolds and the computation of the matrix coe cients in the Galerkin matrix. They satisfy the requirements for a second compression which was introduced to get rid of logarithmic 4 terms in the complexity estimates. We like to point out that the primal multiresolution realizes the order of approximation d, i.e, inf vj2Vj kv vjkL2([0;1]) . 2 jdkvkd; v 2 Hd([0; 1]): The dual multiresolution is generated by compactly supported scaling functions realizing a certain order of approximation e d (d+ e d even). This is not of such importance as in signal analysis, because the dual bases can be avoided in the actual computations since the Galerkin method requires the adjoint and not the dual of the wavelet transform. Nevertheless, there are several situations where the inverse wavelet transform is a very helpful tool, too. With the aid of these scaling functions, we obtain re nable spaces V [0;1] j and e V [0;1] j which contain all polynomials of degree less than d and e d, respectively. The spaces W [0;1] j and f W [0;1] j are de ned uniquely via (4). Biorthogonal wavelet bases [0;1] j = [0;1] j;k : k 2 r[0;1] j ; e [0;1] j = e [0;1] j;k : k 2 r[0;1] j ; generating these complementary spaces, i.e., W [0;1] j := span [0;1] j ; f W [0;1] j := span e [0;1] j ; (5) are not determined uniquely. Each pair of matrices M[0;1] j;1 ;f M[0;1] j;1 2 Rj [0;1] j+1 j jr[0;1] j j satisfying hM[0;1] j;0 M[0;1] j;1 iT hf M[0;1] j;0 f M[0;1] j;1 i = I: de nes wavelets with (5) via [0;1] j = [0;1] j+1M[0;1] j;1 ; e [0;1] j = e [0;1] j+1f M[0;1] j;1 : (6) But, for instance, this straightforward construction does not imply xed and nite masks in the two scale relations of the collections [0;1] j and e [0;1] j . Hence, in order to de ne suitable wavelet bases we utilize the concept of the stable completion [3]. This concept is universal and often employed in the sequel but to avoid confusion we add the su x [0; 1]. De nition 1.1. Let [0;1] j = [0;1] j;k : k 2 r[0;1] j V [0;1] j+1 be a given collection of functions satisfying [0;1] j = [0;1] j+1 M[0;1] j;1 ; M[0;1] j;1 2 Rj [0;1] j+1 j jr[0;1] j j; 5 such that M[0;1] j;0 M[0;1] j;1 is invertible. We de ne the matrix G[0;1] j;0 G[0;1] j;1 with G[0;1] j;0 2 Rj [0;1] j+1 j j [0;1] j j and G[0;1] j;1 2 Rj [0;1] j+1 j jr[0;1] j j as the inverse of M[0;1] j;0 M[0;1] j;1 T , i.e., hM[0;1] j;0 M[0;1] j;1 iT hG[0;1] j;0 G[0;1] j;1 i = I: (7) The collection [0;1] j is called a stable completion of [0;1] j if hM[0;1] j;0 M[0;1] j;1 i l2( [0;1] j+1 ) hG[0;1] j;0 G[0;1] j;1 i l2( [0;1] j+1 ) 1: (8) The stable completion is projected onto W [0;1] j in order to get the desired primal wavelet basis, cf. [8]. In terms of the re nement matrices, the matrix M[0;1] j;1 is de ned by M[0;1] j;1 = hI M[0;1] j;0 f M[0;1] j;0 Ti M[0;1] j;1 =: M[0;1] j;1 M[0;1] j;0 L[0;1] j : (9) One readily veri es that the matrix L[0;1] j 2 Rj [0;1] j j jr[0;1] j j satis es L[0;1] j = f M[0;1] j;0 T M[0;1] j;1 = e [0;1] j ; [0;1] j L2([0;1]): (10) Moreover, one concludes from the identity I = hM[0;1] j;0 M[0;1] j;1 iT hf M[0;1] j;0 f M[0;1] j;1 i = I L[0;1] j 0 I T hM[0;1] j;0 M[0;1] j;1 iT hG[0;1] j;0 G[0;1] j;1 i I 0 L[0;1] j T I the equality hf M[0;1] j;0 f M[0;1] j;1 i = hG[0;1] j;0 +G[0;1] j;1 L[0;1] j T G[0;1] j;1 i ; i.e., f M[0;1] j;1 = G[0;1] j;1 . Remark 1.2. The de nition of M[0;1] j;1 implies [0;1] j = [0;1] j+1M[0;1] j;1 = [0;1] j+1 M[0;1] j;1 [0;1] j+1 M[0;1] j;0 L[0;1] j = [0;1] j [0;1] j L[0;1] j : Consequently, similarly to [31], the wavelets [0;1] j are obtained by updating [0;1] j by linear combinations of the coarse level generators [0;1] j . 6 We abbreviate [0;1] j0 1 := [0;1] j0 ; e [0;1] j0 1 := e [0;1] j0 ; r[0;1] j0 1 = [0;1] j0 1; and de ne the collections [0;1] = [ j j0 1 [0;1] j ; e [0;1] = [ j j0 1 e [0;1] j : Then, according to [8] the following statements hold. The collections [0;1] and e [0;1] de ne biorthogonal Riesz bases in L2([0; 1]). The functions contained in the collections [0;1] and e [0;1] have e d and d vanishing moments, respectively. That is Z 1 0 x j;k(x)dx = 0; = 0; 1; : : : ; e d 1; and likewise for the duals. The functions of the collections [0;1] and e [0;1] have the same regularity as the biorthogonal spline wavelets in L2(R) [33]. Therefore, the norm equivalences are valid in the same ranges, i.e., kvk2Hs([0;1]) X j j0 1 X k2r[0;1] j 2js v; e [0;1] j;k L2([0;1]) 2; s 2 ( e ; ); kvk2Hs([0;1]) X j j0 1 X k2r[0;1] j 2js v; [0;1] j;k L2([0;1]) 2; s 2 ( ; e ); (11) with and e denoting the regularity of the primal and dual wavelets, respectively. Clearly, the goal is to construct a wavelet basis such that only a few boundary wavelets do not coincide with translates and dilates of the Cohen-Daubechies-Feauveau wavelets [4]. For the treatment of boundary integral equations we focus mainly on piecewise constant and linear wavelets, i.e., d = 1 and d = 2. On the level j, we consider the interval [0; 1] subdivided into 2j equidistant subintervals. Then, of course, V [0;1] j is the space generated by 2j and 2j + 1 piecewise constant and linear scaling functions, respectively. The Haar basis and the hierarchical basis on the given partitioning de ne suitable stable completions. Unfortunately, the general situation is not that simple since we have to modify the boundary functions. For the sake of brevity we refer to [8, 21] for the details. The discretization of the hypersingular operator requires globally continuous piecewise linear wavelet bases. According to [14, 21], for their construction, both the 7 scaling functions and the stable completion are required to satisfy certain boundary conditions. The scaling functions [0;1] j and e [0;1] j have to be chosen such that [0;1] j;k (0) = (2j=2; k = 0; 0; k 6= 0; e [0;1] j;k (0) = (2j=2c; k = 0; 0; k 6= 0; (12) c 6= 0, and likewise for x = 1 and k = 2j + 1. This is performed by a suitable change of bases. Additionally, the stable completion [0;1] j is supposed to ful ll zero boundary conditions [0;1] j;k (0) = [0;1] j;k (1) = 0; k 2 r[0;1] j ; (13) as well as the symmetry conditions [0;1] j;k (x) = [0;1] j;3 2j+1 k(1 x); x 2 [0; 1]; k 2 r[0;1] j : (14) For example, the hierarchical basis satis es the latter conditions. 1.3 Wavelets on the unit square In general, it su ces to consider two dimensional wavelets for the treatment of boundary integral equations. Hence, we restrict ourselves to the two dimensional case since the construction keeps simple. For the higher dimensional case we refer to [21]. 1.3.1 Biorthogonal scaling functions The canonical de nition of biorthogonal multiresolutions on the unit square is to take tensor products of the univariate constructions. Hence, the collections of scaling functions are given by j = [0;1] j [0;1] j ; e j = e [0;1] j e [0;1] j ; (15) with the index set j = [0;1] j [0;1] j . Consequently, the associated re nement matrices are M j;0 =M[0;1] j;0 M[0;1] j;0 ; f M j;0 = f M[0;1] j;0 f M[0;1] j;0 : (16) As an immediate consequence of the univariate case, the spaces V j := span j and e V j := span e j are nested and dense in L2( ). Clearly, these spaces are exact of the order d and e d, respectively. We emphasize that the complement spaces W j and f W j are uniquely determined by (4). With this in mind, the remainder of this subsection is dedicated to the construction of biorthogonal wavelet bases j and f j with W j := span j and f W j := span e j . 8 1.3.2 Tensor product wavelets First, we introduce the simplest construction, namely tensor product wavelets j = [0;1] j [0;1] j [ [0;1] j [0;1] j [ [0;1] j [0;1] j : Then, the re nement matrices are de ned via M j;1 = 264M[0;1] j;0 M[0;1] j;1 M[0;1] j;1 M[0;1] j;0 M[0;1] j;1 M[0;1] j;1 375 ; f M j;1 = 264f M[0;1] j;0 f M[0;1] j;1 f M[0;1] j;1 f M[0;1] j;0 f M[0;1] j;1 f M[0;1] j;1 375 : Hence, we di er three types of wavelets on , cf. gures 1 and 2. The rst type is the tensor product [0;1] j;k1 [0;1] j;k2 . The second type is the tensor product of [0;1] j;k1 [0;1] j;k2 . The third type consists of the tensor product of two wavelets [0;1] j;k1 [0;1] j;k2 . We mention that j [0;1] j j jr[0;1] j j implies nearly identical cardinalities of the three types of wavelets. 1.3.3 Simpli ed tensor product wavelets Next, we consider an extension of the tensor product construction. As we will see it replaces the wavelet of the third type by a smoother one. We mention that this simpli es numerical integration, for instance in the Galerkin scheme. The idea is to construct a suitable stable completion on the unit square. Based on the univariate case it can be de ned by the collection j = [0;1] j [0;1] j [ [0;1] j [0;1] j [ [0;1] j [0;1] j : The re nement matrices M j;1, G j;0 and G j;1 are computed by M j;1 = 264M[0;1] j;0 M[0;1] j;1 M[0;1] j;1 M[0;1] j;0 M[0;1] j;1 M[0;1] j;1 375 ; G j;0 = G[0;1] j;0 G[0;1] j;0 ; G j;1 = 264G[0;1] j;0 G[0;1] j;1 G[0;1] j;1 G[0;1] j;0 G[0;1] j;1 G[0;1] j;1 375 : Next, the matrix L j is given by L j = 264 I(j [0;1] j j) L[0;1] j L[0;1] j I(j [0;1] j j) L[0;1] j L[0;1] j 375 : This impliesM j;1 = M j;1 M j;0L j = 264 M[0;1] j;0 M[0;1] j;1 M[0;1] j;1 M[0;1] j;0 M[0;1] j;1 M[0;1] j;1 M[0;1] j;0 M[0;1] j;0 L[0;1] j L[0;1] j 375 : 9 tensor product wavelets 6 x1 x 2 scaling function ? +1 wavelet of type one ? +18 1 +1 18 wavelet of type two ? 18 +1 1 +1 8 wavelet of type three ? + 1 64 18+1 8 1 64 18 +1 1 +1 8 +1 8 1+1 18 1 64 +18 1 8 + 1 64 simpli ed tensor product wavelets 6 x1 x 2 scaling function ? +1 wavelet of type one ? +18 1 +1 18 wavelet of type two ? 18 +1 1 +1 8 wavelet of type three ? + 1 64 1 64 +1 1 1+1 1 64 + 1 64 wavelets optimized with respect to their support 6 x1 x 2 scaling function ? +1 wavelet of type one ? 18 +1 1 +1 8 wavelet of type two ? +1 8 1 +1 18 Figure 1: Interior piecewise constant wavelets with three vanishing moments. The boundary wavelets are not plotted. 10 scaling functions tensor product wavelets wavelet of type one wavelet of type two wavelet of type three simpli ed tensor product wavelets wavelet of type one wavelet of type two wavelet of type three wavelets optimized with respect to their support wavelet of type one wavelet of type two Figure 2: Interior piecewise linear wavelets with four vanishing moments. The boundary wavelets are not plotted. 11 Hence, we di er again three types of wavelets on . The rst and the second type coincide with the tensor product wavelets, see gures 1 and 2. But now the third type consists of the tensor product of the stable completion [0;1] j;k1 [0;1] j;k2 and certain scaling functions j;k0 = [0;1] j;k0 1 [0;1] j;k0 2 of the coarse grid j. In general, the support of this wavelet does not depend on the choice of the stable completion. But choosing a stable completion on [0; 1] with small supports, the product [0;1] j;k1 [0;1] j;k2 2 V j+1 has also small support. Since the additional scaling functions belong to V j , the wavelet is smoother than the corresponding tensor product wavelet. 1.3.4 Wavelets optimized with respect to their supports Last, we consider a more advanced construction which yields wavelets with very small supports. We de ne the wavelet functions via the collections j = [0;1] j [0;1] j [ [0;1] j+1 [0;1] j ; e j = e [0;1] j e [0;1] j [ e [0;1] j+1 e [0;1] j : The re nement matrices M j;1 and f M j;1 are computed by M j;1 = "M[0;1] j;1 M[0;1] j;0 Ij [0;1] j+1 j M[0;1] j;1 # ; f M j;1 = "f M[0;1] j;1 f M[0;1] j;0 Ij [0;1] j+1 j f M[0;1] j;1 # : Thus, we obtain two types of wavelets on , cf. gures 1 and 2. The rst type is the tensor product [0;1] j;k1 [0;1] j;k2 . The second type is the tensor product [0;1] j+1;k1 [0;1] j;k2 . Notice that these wavelets have a very small support in comparison with the previously introduced wavelets, since a scaling function of the ne grid j+1 appears in the rst coordinate. Additionally, the number of wavelets of type two is nearly twice as much as the number of wavelets of type one. Remark 1.3. This construction is highly attractive in higher dimensions since each wavelet coincides only in one coordinate with a wavelet from the interval while in the other coordinates only scaling functions of the levels j and j + 1 appear. 1.4 Wavelets on domains and manifolds In this subsection, we employ a domain decomposition strategy and introduce a family of parametric representations which describe the given manifold. Subsequently, the wavelets on the manifold are de ned via the parametrization. We consider two di erent constructions for wavelets on manifolds. The rst one leads to wavelets which are de ned on each patch individually. The second and more complicated one yields globally continuous wavelets. 12 1.4.1 Parametric representations of manifolds Let denote the unit square, i. e., = [0; 1]2. We subdivide the given manifold 2 R3 into several patches = N [ i=1 i; i = i( ); i = 1; 2; : : : ; N; (17) such that each i : ! i de nes a di eomorphism of onto i. The intersection i\ i0 , i 6= i0, of the patches i and i0 is supposed to be either ; or a common edge or vertex. On the level j, the unit square is subdivided equidistantly j times into 22j squares Cj;k , where k = (k1; k2) with 0 km < 2j. This yields 22jN elements i;j;k := i(Cj;k) i, i = 1; 2; : : : ; N . The construction of globally continuous wavelets requires additionally that the collection of all elements f i;j;kg on a xed level j forms a regular mesh on . Therefore, the parametric representation is subjected the following matching condition. For all x 2 i \ i0 exists a bijective, a ne mapping : ! such that i(s) = ( i0 )(s) = x for s = [s1; s2]T 2 with i(s) = x, cf. gure 3. Unfortunately, this essential supposition restricts the choice of the parametric representation. Figure 3: The parametric representation of the unit sphere is obtained by projecting it onto the cube [ 1; 1]3, which yields six patches (left). On the right hand side one gures out the partition on the level j = 4. The rst fundamental tensor of di erential geometry is given by the matrix Ki(s) 2 R2 2 with Ki(s) := h @ i(s) @sj ; @ i0(s) @sj0 l2(R3)ij;j0=1;2: (18) 13 Since i is supposed to be a di eomorphism, the matrix Ki(s) is symmetric and positive de nite. The canonical inner product in L2( ) is given by (u; v)L2( ) = Z u(x)v(x)d x = N Xi=1 Z u i(s) v i(s) qdet Ki(s) ds: (19) The corresponding Sobolev spaces are indicated by Hs( ). Of course, depending on the global smoothness of the surface, the range of permitted s 2 R is limited to s 2 ( s ; s ). An important role is played by the following modi ed inner product which arises from (19) by omitting the square root of det Ki(s) hu; vi = N Xi=1 (u i; v i)L2( ) = N Xi=1 Z u i(s) v i(s) ds: (20) In L2( ) both inner products de ne equivalent norms hu; ui (u; u)L2( ). However, in general, even on smooth surfaces Ki(s) is not continuous accross the interfaces of the patches, i.e., for i(s) = i0(s0) = x 2 i \ i0 one has Ki(s) 6= Ki0(s0): (21) In the next subsections, the biorthogonal multiresolution on is constructed with respect to the modi ed inner product. Thus, due to (21), the norm equivalence with respect to the canonical Sobolev spaces Hs( ) is limited from below by s = 1=2. 1.4.2 Patchwise smooth wavelet bases If the wavelet basis is not required globally continuous, one may employ wavelet bases de ned on each patch individually. This strategy re ects the canonical one for the piecewise constants. But in the case of piecewise linear functions we obtain multiple de ned knots along the interfaces of intersecting patches. This leads to more degrees of freedom than in the case of global continuity. The primal scaling functions and wavelets are given by i j;k(x) := ( j;k 1 i (x) ; x 2 i; 0; else; i j;k(x) := ( j;k 1 i (x) ; x 2 i; 0; else: Setting i j = i j;k : k 2 j and i j = i j;k : k 2 r j , the collections of scaling functions and wavelets on are de ned by j := i j Ni=1 and j := i j Ni=1. Obviously, the re nement matrices with j = j+1M j;0 and j = j+1M j;1 are obtained by M j;0 = diag M j;0; : : : ;M j;0 | {z } N times ; M j;1 = diag M j;1; : : : ;M j;1 | {z } N times : 14 Clearly, the spaces V j := span j are nested. In addition, we nd V j+1 = V j W j , where W j := span j . Proceeding analogously on the dual side yields a multiresolution on which is biorthogonal with respect to the modi ed inner product (20). Analogously to the univariate case we de ne the collection of wavelets = [ j j0 1 j ; e = [ j j0 1 e j ; with j0 1 := j0 and e j0 1 := e j0. In order to formulate the properties of the wavelets we introduce new function spaces on . For arbitrary s 0 we de ne the Sobolev spaces Hs h ; i( ) as closure of all patchwise C1-functions on with respect to the norm kvkHs h ; i( ) := N Xi=1 kv ikHs( ) : The space L2h ; i( ) indicates as usual the Sobolev space H0 h ; i( ). The Sobolev spaces of negative order, that is H s h ; i( ), are de ned as the duals of Hs h ; i( ) with respect to the modi ed inner product (20). Consequently, H s h ; i( ) is equipped by the norm kvkH s h ; i( ) := sup w2Hs h ; i( ) hv; wi kwkHs h ; i( ) : Proposition 1.4. The collection of wavelets and e form biorthogonal Riesz bases in L2h ; i( ). The primal wavelets have e d vanishing moments in terms of Z s j;k i(s) ds = 0; j j < e d; i = 1; 2; : : : ; N; (22) where = ( 1; 2) denotes a multi index and j j := 1 + 2. Moreover, the norm equivalences kvk2Hs h ; i( ) X j j0 1 X k2r j 22sj hv; e j;ki 2; s 2 ( e ; ); kvk2Hs h ; i( ) X j j0 1 X k2r j 22sj hv; j;ki 2; s 2 ( ; e ); hold with and e corresponding to the Cohen-Daubechies-Feauveau wavelets [4]. We remark that the Sobolev spaces Hs( ) and Hs h ; i( ) are isomorphic in the range s 2 1 2 ; 1 2 , see [14] for details. Therefore, the norm equivalences with respect to 15 the canonical Sobolev spaces are valid in the ranges kvk2Hs( ) X j j0 1 X k2r j 22sj v; e j;k L2( ) 2; s 2 min 1 2 ; e ;min 1 2 ; ; kvk2Hs( ) X j j0 1 X k2r j 22sj v; j;k L2( ) 2; s 2 min 1 2 ; ;min 1 2 ; e : In particular, s = 0 implies the Riesz property of the collection and e , respectively. Remark 1.5. The moment property (22) implies the cancellation property with respect to both inner products. More precisely, let A : Hq h ; i( ) ! H q h ; i( ) denote a boundary integral operator of the order 2q. Then, the estimate hA j0; 0 ; j; i . 2 (j+j0)(e d+1) k dist(supp j; ; supp j0; 0)k2+2e d+2q (23) is satis ed. Analogously, considering a boundary integral operator e A : Hq( ) ! H q( ) of the order 2q we nd ( e A j0; 0; j; )L2( ) . 2 (j+j0)(e d+1) k dist(supp j; ; supp j0; 0)k2+2e d+2q : (24) According to [30], this cancellation property is su cient to compress the system matrices of A and e A, for instance in a Galerkin scheme. 1.4.3 Globally continuous piecewise linear wavelets The construction of globally continuous piecewise linear wavelets on is based on the simpli ed tensor product wavelets. Both, the scaling functions and the stable completion, are required to satisfy the boundary conditions. In order to perform the continuity a gluing technique is utilized along the interfaces of intersecting patches. For the sake of clearness in representation, we introduce rst some further notation since we have to deal with local indices and functions de ned on the parameter domain as well as global indices and functions on the surface . To this end, it is convenient to identify the basis functions with physical grid points on the mesh on the unit square, i.e., we employ the bijective mapping qj : j ! ; qj(k) = 2 jk; in order to rede ne our index sets on the unit square. Then, the boundary conditions (12) and (13) imply j;k @ e j;k @ 0; k 2 j \ ; j;k @ 0; k 2 r j \ : 16 That is, all functions corresponding to the indices k lying in the interior of satisfy zero boundary conditions. Moreover, given any a ne mapping : ! , there holds j = j ; j = j ; e j = e j ; e j = e j : (25) A given point x 2 might have several representations x = i1(s1) = : : : = ir(x)(sr(x)) if x belongs to di erent patches i1 ; : : : ; ir(x). Of course, this occurs only if x lies on a common edge or vertex of these patches. We count the number of preimages to a given point x 2 by the function r : ! N , r(x) := i 2 f1; 2; : : : ; Ng : x 2 i : (26) Clearly, one has r(x) 1, where r(x) = 1 holds for all x lying in the interior of the patches i. Moreover, r(x) = 2 for all x which belong to an edge and are di erent from a vertex. Next, given two points x;y 2 , the function c : ! N de ned by c(x;y) := i 2 f1; 2; : : : ; Ng : x 2 i ^ y 2 i (27) counts the number of patches i containing both points together. Now, on , the index sets are given by physical grid points on the surface j := i(k) : k 2 j ; i 2 f1; 2; : : : ; Ng ; r j := j+1 n j : (28) The gluing along the intersections of the patches is performed as follows. According to [14], the scaling functions j := j; : 2 j and e j := e j; : 2 j are de ned by j; (x) = ( j;k 1 i (x) ; 9(i;k) : i(k) = ^ x 2 i; 0; elsewhere: e j; (x) = ( 1 r( ) e j;k 1 i (x) ; 9(i;k) : i(k) = ^ x 2 i; 0; elsewhere: (29) On the primal side, this de nition re ects the canonical strategy. On the dual side, the strategy is analogously except for normalization, see also gure 4. The normalization factor ensures biorthogonality with respect to the modi ed inner product (20), i.e., h j ; e j i = I. The scaling functions are re nable Riesz bases of the spaces V j := span j and e V j := span e j . The re nement matrices corresponding to these scaling functions j = j+1M j;0; e j = e j+1f M j;0; M j;0; f M j;0 2 Rj j+1 j j j j; 17 XXXXXX 1 2 3 H H Y -1 3 H H Y Figure 4: The primal (left) and the dual (right) generator on a degenerated vertex in the case (d; e d) = (2; 4). are given by M j;0 0; = ( M j;0 k0;k; 9(i;k;k0) : = i(k) ^ 0 = i(k0); 0; elsewhere; (30) f M j;0 0; = ( c( ; 0) r( ) f M j;0 k0;k; 9(i;k;k0) : = i(k) ^ 0 = i(k0); 0; elsewhere: (31) The stable completion j = j; : 2 r j is de ned analogously to the primal scaling functions. In accordance with [14] one has j; (x) = ( j;k 1 i (x) ; 9(i;k) : i(k) = ^ x 2 i; 0; elsewhere: (32) 18 Consequently, the re nement matrix j = j+1 M j;1; M j;1 2 Rj j+1 j jr j j; is determined analogously to M j;0 by M j;1 0; = ( M j;1 k0;k; 9(i;k;k0) : = i(k) ^ 0 = i(k0); 0; elsewhere: (33) The dual wavelets e j := e j; : 2 r j are obtained by their re nement relation. There holds e j = e j+1f M j;1; f M j;1 2 Rj j+1 j jr j j; with f M j;1 given by f M j;1 0; = ( c( ; 0) r( ) f M j;1 k0;k; 9(i;k;k0) : = i(k) ^ 0 = i(k0); 0; elsewhere; (34) cf. [14, 21]. A closed expression of the matrix L j := f M j;0 T M j;1 2 Rj j j jr j j is crucial for the implementation of the discrete wavelet transform. In fact, the next theorem con rms the existence of an explicit formula [21]. Theorem 1.6. The entries of the matrix L j = f M j;0 T M j;1 are given by L j 0; = ( c( ; 0) r( 0) L j k0;k; 9(i;k;k0) : = i(k) ^ 0 = i(k0); 0; elsewhere: (35) As an immediate consequence of this theorem, a black box algorithm for the computation of the discrete wavelet transform is available. Of course, the de nitions of the re nement matrices seem to be very technical. However, as the algorithms 1 and 2 con rm, the implementation of the discrete wavelet transform is rather canonical. The collection of wavelets are given by = [ j j0 1 j ; e = [ j j0 1 e j ; with j0 1 := j0 and e j0 1 := e j0 . For arbitrarily s 0 let the Sobolev spaces Hs h ; i( ) be the closure of all globally continuous, patchwise C1-functions on with respect to the norm kvkHs h ; i( ) := N Xi=1 kv ikHs( ) : 19 The space L2h ; i( ) indicates as usual the Sobolev space H0 h ; i( ). The Sobolev space H s h ; i( ), s < 0, indicates the dual of Hs h ; i( ) with respect to the modi ed inner product (20), i.e., kvkH s h ; i( ) := sup w2Hs h ; i( ) hv; wi kwkHs h ; i( ) : Proposition 1.7. The collection of wavelets and e form biorthogonal Riesz bases in L2h ; i( ). The primal wavelets satisfy the cancellation property (23). Furthermore, there holds the norm equivalence kvk2Hs h ; i( ) X j j0 1 X 2r j 22sj hv; e j; i 2; s 2 ( e ; ); kvk2Hs h ; i( ) X j j0 1 X 2r j 22sj hv; j; i 2; s 2 ( ; e ); with and e denoting the regularity of the Cohen-Daubechies-Feauveau wavelets. Unfortunately, the cancellation property is not satis ed with respect to the canonical inner product. Clearly, (24) is valid if both wavelets, j; and j0; 0, are supported on a single patch. We mention that for wavelets supported on several patches this estimate remains true only if the fundamental tensor of di erential geometry Ki is continuous across the interfaces of intersection patches. According to [14], the Sobolev spaces Hs( ) and Hs h ; i( ) are isomorphic in the range s 2 1 2 ;min 32 ; s , This yields the norm equivalences kvk2Hs( ) X j j0 1 X 2r j 22sj v; e j; L2( ) 2; s 2 min 1 2 ; e ;min 32 ; ; kvk2Hs( ) X j j0 1 X 2r j 22sj v; j; L2( ) 2; s 2 min 1 2 ; ;min 32 ; e : In particular, s = 0 implies the Riesz property of the collections and e , respectively. 20 Algorithm 1 This algorithm computes the two scale decomposition e j+1a(j+1) = e j a(j) + e j b(j), where a(j) = a(j) 2 j and b(j) = b(j) 2r j . initialization: a(j) := b(j) := 0 for i = 1 to N do begin for all k 2 j do begin C: compute coe cients of e j for all k0 2 j+1 do begin a(j) i(k) = a(j) i(k) + M j;0 k0;ka(j+1) i(k0)=r i(k) end end for all k 2 r j do begin C: compute coe cients of e j for all k0 2 j+1 do begin b(j) i(k) = b(j) i(k) + M j;1 k0;ka(j+1) i(k0)=r i(k) end end end 21 Algorithm 2 This algorithm computes the two scale decomposition j+1a(j+1) = j a(j) + j b(j), where a(j) = a(j) 2 j and b(j) = b(j) 2r j . initialization: a(j) := b(j) := 0 for i = 1 to N do begin for all k 2 j do begin C: compute coe cients of j for all k0 2 j+1 do begin a(j) i(k) = a(j) i(k) + M j;0 k0;ka(j+1) i(k0)=r i(k0) end end for all k 2 r j do begin C: compute coe cients of j for all k0 2 j+1 do begin b(j) i(k) = b(j) i(k) + M j;1 k0;ka(j+1) i(k0)=r i(k0) end end end for i = 1 to N do begin for all k 2 r j do begin C: add scaling functions to j for all k0 2 j do begin b(j) i(k) = b(j) i(k) L j k0;ka(j) i(k0)=r i(k0) end end end 22 2 The wavelet Galerkin scheme This section presents a fully discrete wavelet Galerkin scheme for boundary integral equations. In the rst subsection we discretize a given boundary integral equation. Then, in subsection 2.2 we introduce the a-priori matrix compression which reduces the relevant matrix coe cients to an asymptotically linear number. In subsections 2.3 and 2.4 we point out the computation of the compressed matrix. Next, in subsection 2.5 we introduce an a-posteriori compression which reduces again the number of matrix coe cients. The last subsection is dedicated to the preconditioning of system matrices which arise from boundary integral operators of nonzero order. For the sake of simplicity, we skip the su ces of the spaces and their bases. The collection J := SJ 1 j j0 1 with a capital J denotes the wavelet basis of VJ . Moreover NJ := j J j indicates the number of unknowns on the level J . 2.1 Discretization Generally written, a boundary integral equation for the unknown density 2 Hq( ) is given by A = f on : (36) Hereby, A : Hq( ) ! H q( ) denotes a boundary integral operator of the order 2q and f 2 H q( ) indicates the right hand side. Consequently, the variational formulation of (36) reads seek 2 Hq( ) : (A ; )L2( ) = (f; )L2( ) 8 2 Hq( ): (37) It is well known, that the variational formulation (37) is equivalent to the boundary integral equation (36), see e.g. [18, 25] for details. For the Galerkin scheme we replace the energy space Hq( ) in the variational formulation (37) by the nite dimensional spaces VJ introduced in the previous section. Then, we arrive at the problem seek J 2 VJ : (A J ; J)L2( ) = (f; J)L2( ) 8 J 2 VJ : Equivalently, due to the nite dimension of VJ , the ansatz J = J J together with A J := A J ; J L2( ); f J := f; J L2( ); yields the wavelet Galerkin scheme A J J = f J : (38) The system matrix A J is quasi-sparse and might be compressed to O(NJ) nonzero matrix entries if the wavelets have a su ciently large number of vanishing moments. The remainder of this paper is devoted to the e cient computation of the system matrix. 23 Remark 2.1. Replacing the wavelet basis J by the single-scale basis J yields the traditional single-scale Galerkin scheme A J J = f J ; where A J := A J ; J L2( ), f J := f; J L2( ) and J = J J . This scheme is related to the wavelet Galerkin scheme by A J = TJA JTTJ ; J = T T J J ; f J = TJf J ; where TJ denotes the wavelet transform. The system matrixA J is densely populated. Therefore, the costs of solving a given boundary integral equation traditionally in the single-scale basis is at least O(N2 J). 2.2 A-priori compression The discretization of a boundary integral operator A : Hq( )! H q( ) by wavelets with a su ciently large number of vanishing moments (22) or a corresponding cancellation property (23) yields quasi-sparse matrices. In a rst compression step all matrix entries, for which the distances of the supports of the corresponding ansatz and test functions are bigger than a level depending cut-o parameter Bj;j0, are set to zero. In the second compression step also some of those matrix entries are neglected, for which the corresponding ansatz and test functions have overlapping supports. First, we introduce the abbreviation j;k := conv hull(supp j;k); j;k := sing supp j;k: Note that j;k denotes the convex hull to the support of j;k while j;k denotes the so-called singular support of j;k, i.e., those points where j;k is not smooth. The compressed system matrix A J corresponding to the boundary integral operator A is de ned by [A J ](j;k);(j0;k0) := >>><>>>: 0; dist j;k; j0;k0 > Bj;j0; j; j 0 j0; 0; dist j;k; j0;k0 > B0 j;j0; j 0 > j; 0; dist j;k; j0;k0 > B0 j;j0; j > j 0; A j0;k0; j;k L2( ); otherwise: (39) Herein, choosing a; a0 > 1; d < ; 0 < e d+ 2q; (40) 24 the cut-o parameters Bj;j0 and B0 j;j0 are set as follows Bj;j0 = a maxn2 minfj;j0g; 2 2J( q) (j+j0)( +e d) 2( e d+q) o; B0 j;j0 = a0maxn2 maxfj;j0g; 2 2J( 0 q) (j+j0) 0 maxfj;j0g e d e d+2q o: (41) The resulting structure of the compressed matrix is guratively called nger structure, cf. gure 5. It is shown in [30] that this compression strategy does not compromise the stability and accuracy of the underlying Galerkin scheme. Theorem 2.2. Let the system matrix A J be compressed in accordance with (39), (40) and (41). Then, the wavelet Galerkin scheme is stable and the error estimate k JkH2q d( ) . 2 2J(d q) Hd( ) (42) holds, where 2 Hd( ) denotes the exact solution of the given boundary integral equation A = f and J = J J is the numerically computed solution, i.e., e A J J = f . Consequently, we obtain the optimal order of convergence of the Galerkin scheme. 0 100 200 300 400 500 600 700 800 90